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The top quark antiquark production system in the dilepton decay channel is described by a set 
of equations which is nonlinear in the unknown neutrino momenta. Its most precise and least time 
consuming solution is of major importance for measurements of top quark properties like the top 
quark mass and tt spin correlations. The initial system of equations can be transformed into two 
polynomial equations with two unknowns by means of elementary algebraic operations. These two 
polynomials of multidegree two can be reduced to one univariate polynomial of degree four by means 
of resultants. The obtained quartic equation is solved analytically. 



PACS numbers: PACS29.85.+C 
I. INTRODUCTION 



In 1992, Dalitz and Goldstein published a numerical 
method based on geometrical considerations to solve the 
system of equations describingthe kinematics of the tt 
decay in the dilepton channel [l[ 0] • In 2004 an approx- 
imation of the system of equations - assuming that the 
transverse momentum of the tt system can be neglected 
- has been solved analytically Q by means of computer 
algebra software such as Meanwhile the transverse 
momentum constraint has been omitted while the solu- 
tion is still derived by means of computer algebra and 
its accuracy does not reach real precision [Fj. In 2005, 
the system of equations could be solved algebraically to 
real precision free of any singularity [||. The analyti- 
cal solution introduced here is based on a new Ansatz 
which minimises the amount of intermediate steps to de- 
rive the solution. This approach makes the need of com- 
puter algebra superfluous. In addition it provides more 
transparency and control over singularities which are in- 
trinsic to the analytical solution. Further the accuracy 
achieved is - as already in the algebraic approach [|| - 
of real precision. Important improvements in terms of 
robustness, code volume and time consumption with re- 
spect to the algebraic approach make this method more 
convenient for applications in practice. Other solution 
methods can compare their performance to the reference 
method described here. It should be mentioned that dif- 
ferent approachs leading to analytical solutions, without 
giving a complete algebraic derivation and without rig- 
orous discussion of reducible and irreducible singularities 
exist in the literature 0] @|. 

In the next section the system of tt dilepton equations 
is introduced, followed by the derivation of the analytical 
solution including a rigorous discussion of the reducible 
and irreducible singularities of the analytical solution. 
Subsequently, the performance of the method is elabo- 
rated. 



II. tt DILEPTON KINEMATICS 

The system of equations describing the kinematics of tt 
dilepton events can be expressed by the two linear and 
six non linear equations 

$x = Pu x + Pv x , 
$y = Pv v +Pu v , 

El = Pl x +pl y + 
El = V% x +Pl v + pl m , 
m 2 w + = (E i+ + E v ) 2 - (p e + +p„ m ) 2 , 

~(Pe+ +Pv v ) 2 - (Pe+ +PuJ 2 , 
m 2 w- = {Et- + E,f - ( Pe - +p p J 2 , (1) 

~(Pe- +Pp v ) 2 - (Pe- +PpJ 2 : 
m 2 = {E b + E e + + E v f - {p K + Pe + +p„J 2 , 

-fay +P £ + +PVy) 2 - fa. +P t + +PuJ 2 , 

ml = {E h + E t - + E v f - (jr bx p t - +p p J 2 , 

~(Pby +P t - +Pvyf - (Pb. +P(- +P*J 2 - 

The z-axis is here assumed to be parallel orientated to the 
beam axis while the x- and y-coordinates span the trans- 
verse plane. The first two equations relate the projection 
of the missing transverse energy onto one of the trans- 
verse axes (x or y) to the sum of the neutrino and an- 
tineutrino momentum components belonging to the same 
projection. The next two equations relate the energy 
of the neutrino and antineutrino, which are assumed to 
be massless in good approximation, with their momenta. 
Finally four non linear equations describe the W boson 
and top quark (antiquark) mass constraints by relating 
the invariant masses to the energy and momenta of their 
decay particles via relativistic 4- vector arithmetics. 

III. ANALYTICAL SOLUTION 
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The system of equations (|T|) can be subdivided in two 
entangled sets of equations. One set of equations, de- 
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scribing the t — > bW + — > M + i/f parton branch of the 
event, depends on p Vz while the other pair of equations, 
describing the i — > bW~~ — > M~ parton branch of the 
event, depends on p Pz . 

The equation describing the invariance of the W boson 
mass can be expressed in the following way 



(E l+ +E v ) 2 -{p7 + +pl) 2 
Ej + + 2Eg+ E v + E 2 . 



*+ -2-pi+fv-fv 



2E t + E v — 2pi+pl 



which can be rewritten as 



E v = 



b w+ 



2E, 



(+ 



(2) 



(3) 



The equation describing the invariance of the top quark 
mass can be transformed in the same way leading to 



E v = - 



-2E b E e ++2p b p ( '+ +2(p b +pe+)p~l 



2(E b + E i+ ) 



(4) 



where additional terms emerge due to the fact that quan- 
tities which depended in equation {3]) only on the lepton 
depend now also on the b quark. Next the unknown E v 
can be eliminated by subtracting equation ((4]) from ([3]), 
leading to an equation of the form 







at + a 2 p Vx + a-zpv + a 4 p Uz 



(5) 



where the coefficients a are constants given in the ap- 
pendix. This equation is linear in the three neutrino 
momentum components. Since the unknown p„ a does 
only appear in the top quark parton branch it is manda- 
tory to eliminate this variable with a linear independent 
equation of the top quark parton branch to obtain finally 
together with the equations of the antitop quark branch 
two equations of the two unknowns p Ux and p v . 

To eliminate the unknown p Vr it is straight forward 
to use equation (for convenience multiplied by the 
denominator 2Eg+). The neutrino energy E v can be ex- 
pressed in terms of the three neutrino momenta compo- 
nents in substituting it with the third equation of (fTJ). To 
obtain a polynomial equation the squared of this equa- 
tion is being considered in the following. Its longitudinal 
neutrino momenta can be substituted by equation l|5|). 
solved to p Vz . The resulting equation of the form 

= C22 + C2lPv„+CnPvy+C2lpl a +CioPv x Pvy+CooPl 1/ (6) 

is a multivariate polynomial of multidegree two which 
depends only on the transverse neutrino momenta p Ux 
and p v . The coefficients are again constants which can 
be expressed in terms of the former derived constants a 
and are given in the appendix. 

In the same way can be proceeded for the equations 
describing the antitop quark parton branch. The equiv- 
alent of equation ([5]) reads 



= &i + b 2 pu x + hPi> v + hp*. 



(7) 



and the counter part of polynomial ([6]) can be written as 

= 42 + 4iPp,+4iPp H +4iPL+4oP^Pp ! ,+ d ooPp B ( 8 ) 

The two equations linear in the three (anti-)neutrino 
momenta and (J7J build the minimal Ansatz used here. 
In contrast the Ansatz made in 0] is based on two 
equations linear in the four unknowns p Pjc , pp y , p„ z , p v% . 

To reduce equations ^ and §8§ to two polynomial 
equations of two unknowns the transverse antincutrino 
momenta of equation |[8]) can be expressed by the trans- 
verse neutrino momenta with help of the missing trans- 
verse energy relations of the system of equations ([1}. 
Since these relations are linear in the neutrino and an- 
tineutrino momenta the substitution leads again to a 
polynomial of the form 

= d 2 2+d2ip V:I: +d n p Uy +d2ipl ai +dioPu x Pu y +doopl v (9) 

with multidegree two whose coefficients are given in the 
appendix. To solve these two polynomials without loss 
of generality to p^ the resultant with respect to the neu- 
trino momentum p v is computed as follows. The coef- 
ficients and monomials of the two polynomials ^ and 
([§]) are rewritten in such a way that they are ordered in 
powers of p v like 



coP„ y 
d opL 



C X P Vy 

■ d\p Vv 



C 2 , 

- d 2 



(10) 
(11) 



where c and d are polynomials of the remaining unknowns 
p„ x , p v and the coefficients c rn , d n are univariate poly- 
nomials of p Ux only. The resultant can then be obtained 
by computing the determinant of the Sylvester matrix 



Res(p„ ) = Det 




(12) 
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FIG. 1: Number of solutions per event for particles before any 
radiation, t quark and W boson masses are assumed to be 
known exactly. 
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FIG. 2: Distribution of the expression c\do — codi which ap- 
pears in the denominator in the solution of p Vy . Since the 
distribution is symmetric around zero the module of the ex- 
pression is plotted. As can be seen the values assumed by 
the expression are far away from zero which would cause a 
singularity in the solution. 



As shown in Fig. [2] the coefficient in the denominator of 
equation (|14[) does not acquire values which are even close 
to the singularity at zero. Thus it is ensured that the neu- 
trino momenta p Uy and p„ y can be computed accurately 
over the whole phase space of possible solutions. 

Finally the longitudinal (anti-)neutrino momenta p Vz 
and p Pz can be easily obtained by the linear equations 
and ([7]) assuming that the coefficients 0,4 and 64 
are different from zero since they appear as a prod- 
uct together with the longitudinal (anti-) neutrino mo- 
menta themselves. The distributions of the coefficients 
are shown in Fig. [3] The fraction of solutions close to the 
singularity - irreducible in the analytical solution - is be- 
low the per mill level and may be neglected for practical 
purposes. From a theoretical point of view this singular- 
ity can be circumvented in solving the neutrino momenta 
p Uz and pv z analytically with the equations (2) and (3) 
of the algebraic approach Q which does not contain any 
singularity. It has been verified that the longitudinal 
(anti-)neutrino momentum does not typically vanish to- 
gether with the coefficient 0,4 (64) simultaneously which 



which is equated to zero. The omitted elements of the 
matrix are identical to zero. The resultant is a univariate 
polynomial of the form 



x10" 



h oP^ + hip v + h 2 p„ + h z p Vj: + h 4 



(13) 



which contains the remaining unknown p„ x . It is of de- 
gree four and can be solved analytically. The coefficients 
h are given in the appendix. This result shows that there 
is at most a four fold ambiguity. In Fig.[T]the distribution 
of the number of solutions per event is plotted. Here it 
has been assumed that the 4- vectors of the particles and 
the top quark and W boson masses which enter into the 
system of equations are known exactly. Under these con- 
ditions there are two solutions in about 80% of cases and 
four solutions else. In the next section it will be inves- 
tigated how this distribution changes under more realis- 
tic conditions when the assumption of exactness between 
particles and reconstructed objects is not valid anymore. 
Once the solution of a neutrino momentum p Ux has been 
found the other neutrino and antineutrino momentum 
components have to be determined. The antineutrino 
momentum pp a can be immediately obtained by the lin- 
ear transverse missing energy relation of the initial sys- 
tem of equations ([1]) . To derive the neutrino momentum 
p Vy equation (fTTTj) is multiplied by do and equation (fill) 
is multiplied by cq so that their difference yields a linear 
equation in the neutrino momentum p Uy which can then 
be isolated as 



cpd 2 - c 2 d 
cirf - c di 



(14) 



Again the antineutrino momentum p Py can be imme- 
diately obtained by the corresponding linear transverse 
missing energy relation of the initial system of equations. 



S 0.008 

w 

_o 

I 0.006 
0.004 
0.002 



Entries 905590 
Mean -0.1148 
RMS 2.934 




-2 



X10" 



S 0.008 

W 

o 

I 0.006 
0.004 
0.002 



Entries 905590 
Mean 0.008098 
RMS 2.846 




-2 



FIG. 3: Distributions of the coefficients 04 (top) and 64 (bot- 
tom). The coefficients are flat distributed over the whole 
phase space including the value zero where an irreducible sin- 
gularity of their reciprocal resides. 



would cause the singularity to disappear. jv so1 >o <N soi > R- M S(iV s >|) 



IV. PERFORMANCE OF THE METHOD 

The performance studies discussed here are assuming 
Tevatron proton antiproton collider settings with a centre 
of mass energy of 1.96 TeV which has been set up in the 
Monte Carlo event generator PYTHIA 6.220 0. Cross- 
checks at a centre of mass energy of 14 TeV assuming the 
LHC proton collider environment confirm the indepen- 
dence of the method of particular collider settings. 

The quartic equation in y Vx (|13|) is typically flat around 
the neutrino momenta of interest where solutions are to 
be expected. Fig. [4]shows the function for a given event. 
Only in the bottom plot where the function is zoomed 
out the roots can be recognized. Computationally the 
solutions are robust. The generated (anti-)neutrino mo- 
menta coincide with one of the solutions to real precision 
assuming that the 4- vectors of the two leptons and the b, 
b quarks and the masses of the (anti-)top quarks and W 
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FIG. 4: A typical quartic equation whose real roots in p„ x 
are solutions of the initial system of equations describing the 
ti dilepton kinematics. The bottom plot is zoomed around 
the interesting p Vm range of the abscissa where the analytical 
solution becomes singular. 
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t, W pole mass assumed 
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t, W pole mass assumed 
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0.59 
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1.35 
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0.79 
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(parton matched) 
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3.19 


1.47 


jets + leptons smeared) 
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jets + leptons smeared), 


0.0072 


7.96 


4.72 


reconstructed objects 








100 x resolution smeared 









TABLE I: Number of solutions, fractions and statistical quan- 
tities for events which have been solved (JV so j > 0). The left 
column shows the fraction of events having exactly two solu- 
tions. In the centre the average number of solutions per solved 
event is given. To the right the RMS of this distribution is 
shown. 



bosons which are entered into the quartic equation are 
known exactly. The fraction of events where no solution 
can be found or no solution coincides with the generated 
(anti-)neutrino momenta to real precision is below the 
per mill level. If the W boson mass is generated off- 
shell while its pole mass is assumed in the solution the 
efficiency drops to 89%. Relaxing the same assumption 
for the top quark mass results into a further decrease of 
efficiency to 84%. Beyond, an infrared-safe cone algo- 
rithm flOj with cone size R = 0.5 in the space spanned 
by pseudorapidity and azimuthal angle has been applied 
to the hadronic final state particles. Two reconstructed 
jets, two leptons and missing transverse energy are re- 
quired for an event to be selected. The jets are accepted 
as 6-tagged if they coincide within AR < 0.5 with the 
generated b quarks. The solution efficiency drops to 71% 
and can be re-established at 81% in solving both b quark 
jet permutations. Smearing the leptons and jets with the 
energy resolution of the DO detector [ll| decreases the ef- 
ficiency to 75%. In practice, a given event can be solved 
repeatedly, with the energy of the particles and objects 
smeared randomly within the detector resolution once 
each iteration to improve the solution efficiency. These 
observations are consistent with the findings of the alge- 
braic approach Q. This confirms on one hand the reli- 
ability of the algebraic approach and rises on the other 
hand the question what numerical methods with a supe- 
rior solution efficiency are actually solving. 

Considering only events which could be solved it is 
important to investigate the number of solutions in de- 
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pendence of the experimental settings since this number 
is directly proportional to the ambiguities of the solved 
and reconstructed events which in turn determines the 
significance of the solutions and any observable making 
use of it. In Tab. [I] the fraction of solved events having 
exactly two solutions, the average number of solutions 
and its RMS is given for different experimental settings. 
The first four lines describe the evolution of these quan- 
tities derived from the particle final state. Relaxing the 
amount of assumption about the top quark and W boson 
masses increases the fraction of solved events with exactly 
two solutions while the average number of solutions and 
its RMS decrease slightly. Allowing both b quark jet per- 
mutations - assuming that the charge of the quarks can 
not be determined with adequate certainty - the frac- 
tion of events having exactly two solutions drops con- 
siderably in favour of a higher solution multiplicity with 
larger RMS. The table items below show the number of 
solutions for reconstructed objects, first for right, wrong 
and both b quark jet permutations then energy resolu- 
tion smearing is applied to the reconstructed objects and 
finally 100 solution iterations have been accomplished to 
take into account the uncertainty in the measured en- 
ergy of the reconstructed objects. The general tendency 
is that the fraction of solved events with exactly two so- 
lutions decreases with less accurate knowledge about the 
particles and objects while the solution multiplicity and 
its RMS does increase. 



V. CONCLUSIONS 

The analytical solution of the system of equations de- 
scribing the ti dilepton kinematics has been presented. 
The Ansatz of formulating two equations linear in 
the three neutrino and antineutrino momentum compo- 
nents leads after substitution of the longitudinal (anti- 
neutrino momenta to two multivariate polynomials of 
two unknowns with multidegree two. It turns out that 
each of these two polynomials has a singularity which can 
be removed. In contrast there are two irreducible singu- 
larities in the linear equations described above which can 
be circumvented in exploiting the analytical Ansatz of 
the algebraic approach Q to determine the longitudinal 
(anti-) neutrino momenta. The two multivariate polyno- 
mials can be reduced to a univariate polynomial of de- 
gree four by means of resultants. The obtained quartic 
equation is solved analytically. The solution could be de- 
rived without any use of computer algebra software. The 
fraction of events without any solution or with no so- 
lution matching the generated (anti-)neutrino momenta 
with real precision are below the per mill level assuming 
that the particle momenta and masses inserted into the 
analytical solution are known exactly. Consistent with 
the observations of the algebraic approach Q little devi- 
ations of the inserted particle momenta and masses from 
their true values drop the solution efficiency and purity 
considerably. At the same time the solution multiplicity 



increases. This raises the question what more efficient 
numerical methods are actually solving. General solu- 
tion methods can compare their performance with the 
analytical solution described here. 
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APPENDIX 
Polynomial coefficients 

The coefficients of equation §5§ are given by 

oi = (E b + E t +)(rriw - m 2 l+ ) 
—E e +(m 2 — ml — m 2 +) 
+2E b E 2 e+ -2E i+ f bP i+, 

a 2 = 2(E b p e + - E e +p bx ), 

a 3 = 2(E b p e + - E e +p by ), 

a 4 = 2{E b p t + - E t +p bz ) 

where it is important that the coefficient 04 does not van- 
ish since equation © has to be divided by it to isolate the 
unknown p^. As explained in section HTT1 this irreducible 
singularity can be circumvented in solving for p Vi with 
the analytical Ansatz made in the algebraic approach [|| . 

The equivalent equation of the antitop quark parton 
branch is 

= bi + b 2 p^ + b 3 p Py + bip 9z (15) 
with the coefficients 

bi = (E b + Et-)(m,w - mf_) 

771 / 2 2 2 \ 
-E t -(m t -m- b - m e -) 

+2E b El -2E t -f0>T-, 
b 2 = 2(E b p t - - Ei-p b J, 
b 3 = 2{E b p l - -E e -p by ), 
b 4 = 2{E b p l - - E e -p b J. 

Again there is a singularity in case of vanishing coefficient 
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&4- The coefficients of equation ([6]) are given by 

C22 = {m 2 v+ - m 2 + ) 2 - A(E t + 2 - p 2 +)(ai/a 4 ) 2 

-4(m^ + - mj + )p e +a 1 /a 4 , 
c 2 i = 4(m^+ — m|+)(p^+ -p i +a 2 /a i ) 

-8(E 2 + -p 2 +)aia 2 /a\ - 8p e +p e +ai/a 4 , 

C 20 = -4(£?+ -p^)_ 4 (£;2 + -J^.)(02/04) 2 

-8p e +p e +a 2 /a 4 , 
en = 4(m^+ - m 2 +){p e + -p gas/a*) 

-8(E$+ -p 2 e +)aia3/al - 8p e +p e +ai/a 4l 
cio = ~8{E 2 + -p 2 +)a 2 a 3 /al + 8p e +p t + 

-8p g +p e +a 3 /a 4 - 8p e +p £ +a 2 /a 4 , 
c 00 = -i(E 2 + -p 2 t )-4(E 2 + -p 2 t )(a 3 /a 4 ) 2 

-8p e +p e +a 3 /a 4 . 

Similar the coefficients d! of the antitop quark branch 
depend on the coefficients 6 in the following way 

d 22 = (m 2 w - -m 2 .) 2 -4(E 2 . - p^h/bif 

-4(m^- - m 2 -)p e -b 1 /b 4 , 
d' 21 = 4(m^_ - mj-)(p e - -p t - 62/64) 

-8(^ 2 _ - p 2 -)b x b 2 /b\ - 8 Vl -p r b 1 jb 4 , 

4o = -4(4--Pl-)- 4 (^ 2 --Pl-)(W 2 

~8p e -P e7 b 2 /b 4 , 
d' n = 4(m^_ -m 2 -)(p e - -^-63/64) 

-8{E 2 . - p 2 )hb 3 /bl - 8p e -Pt- 61/64, 
d! w = -8(El-p 2 7 )b 2 b 3 /b 2 4 + 8p e -p e - 

-8p i -p e -b 3 /b 4 - 8p e -p e7 b 2 /b 4 , 
d' 00 = -A{El - p 2 ) - A{El - P 2 -)(b 3 /b 4 ) 2 

-8p l -p i -b 3 /b 4 . 

The remaining unknowns in these equations - which are 
the transverse antineutrino momenta - are substituted by 
the missing transverse energy relations of the system of 
equations ([T]) to obtain finally the set of equations 

d 22 = d' 22 + l£ x 2 d 2Q + lfiy 2 d' 00 + J/i x Jfiyd! X Q 

+Ifixd 2 i + Ifjyd l lll 
d 2 \ — — d 21 — 2J/} x d 20 — lp y d 10 , 
d 2 o = d' 20 , 

dn = —d n — 2I/J y d' 0Q — tfxd' w 
dw = d' w , 
d<x) = d' 00 , 

which depends merely on the transverse neutrino mo- 
menta p v% and p Vy . 



The resultant expressed in terms of the multivariate 
polynomials Cjk and d mn are given by 

h 4 = cl dl 2 + Cnd 22 (ciid 00 - coo^ii) 

+cooc 22 (dj 1 -2dood 22 ) + c 22 d o(c 22 d o- cudn) , 
h 3 = c 00 d 2 i(2c od 22 -ciidii) + c 00 du(2c 22 d w + c 21 dii) 

+c 22 d 00 (2c 2 idoo- ciidro) - c 00 d 22 (ciid w + c w dn) 

-2c od o(c 22 d 2 i + c 2 id 22 ) - doorfii(ciiC2i+cioc 2 2) 

+cndoo(ciid 2 i + 2ci d 22 ), 
h 2 = (? m (2d 22 d 2a + d 2 21 ) - c 00 d 21 (ciid w + c 1Q dn) 

+cnd 2 o(cii(ioo - coodn) + c oc?io(c22dio - Ci rf 22 ) 

+c orfii(2c 2 irfio + c 20 dn) + (2c 22 c 20 + c 21 )do 

-2c odoo(c 2 2d 2 o + c 2 id 2 i + c 20 d 22 ) 

+ciorfoo(2ciid 2 i + c w d 22 ) - d Q0 d w (ciic 2 i + c w c 22 ) 

-doodu(cnc 20 + C10C21), 
hi = cwd 2 i{2c m d 2Q -c w d w ) - c m d 2 Q{cud w + c w dn) 

+ C00^10 (C21 dlO + 2c 2 odll ) - 2c 0<ioO (C21 ^20 + C20^21 ) 

+ciodoo (2cn d 20 + c w d 2 i ) - c 20 d o (2c 2 id 00 - Cio<2ii) 
-doodw(cnc 20 + C10C21), 

^0 = Cqq^O + c 10^2o(ciorfoO — Cooliio) 

+C2orfio(coo^io — ciodoo) + C2odoo(c20^oo — 2000^20)- 

To avoid singularities which arise in the case of van- 
ishing factors a 4 or 64 the coefficients Cjk, d mn of the 
polynomials (J6j> and (|8j have been multiplied with the 
least common multiple of the denominators which are a 2 
and 64 respectively. These factors are constant for a given 
event and thus do not alter the position of the real roots 
which correspond to the neutrino momenta p Vx . 

Quartic equation 

The quartic equation can be solved analytically in re- 
ducing it to a cubic equation. There are several ways to 
achieve this. Here the method of Ferrari [I2J - who was 
the first to develop an algebraic technique for solving the 
general quartic equation - is being used. 

First the leading coefficient ho of the quartic polyno- 
mial (|13p is normalized to one (in the case the leading 
coefficient vanishes the problem is already reduced to a 
cubic equation). If the constant h 4 vanishes the quar- 
tic polynomial can be factorized into p Vx times a cubic 
equation. In this case one root namely p v=c = is already 
known. The substitution p v% — p' v — h±/A leads to the 
simplified equation 

= p'l + kip' 2 x + k 2 p' Vx + k 3 

with the coefficients 

k-i = h 2 -3/if/8, 

k 2 = h 3 + hl/8-hih 2 /2, 

k 3 = hi- 3/^/256 + hlh 2 /16-h 1 h 3 /4:. 
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If the coefficient k 3 vanishes again the equation can be 
factorized into p' v times a cubic polynomial. If the coef- 
ficient k 2 vanishes the quartic polynomial in p' v can be 
expressed as a quadratic equation in p' 2 . In the general 
case where all three coefficients fci , k 2 and k 3 are different 
from zero the quartic polynomial can be factorized into 
the product of two quadratic polynomials as follows 



'4,; 12 

Pu x + k lP Vlc 

(p'Z+tip' v +t 2 )(p'l - 



f- fap'vn + k 3 = 

hp' + h/t 2 ). 



(16) 



Once the new coefficients t\ and t 2 have been determined 
the quadratic polynomials can be easily solved. Compar- 
ison of the coefficients yields 



sign changes. A cubic equation with real coefficients has 
always either one or three real roots. In the case of two 
or zero negative roots there must conclusively be at least 
one positive root. Once this root has been determined, t\ 
can be inserted into equation (|18[) above to determine t 2 
and subsequently the quadratic polynomials (I16p of the 
quartic equation. 

Cubic equation 

There are several ways to solve the cubic equation [3] ■ 
Here the approach of [l5j has been adopted. The cubic 
equation 

= z 3 + siz 2 + s 2 z + s 3 



and 



h = k 3 /t 2 + t 2 -t\ 



h{h/t 2 -t 2 ). 



is assumed to have real coefficients. First the two vari- 
ables 

s\ - 3s 2 



It is ensured that t 2 which appears in the denominator 
does not vanish since the coefficient k 3 has been assumed 
to be different from zero and k\ , k 2 are finite. Eliminating 
t 2 in the two nonlinear equations above yields to a cubic 
equation in t\. To achieve this the two equations above 
are rewritten in the following form 

k 3 /t 2 +t 2 = ki +t\ + k 2 /t u 
k 3 /t 2 - t 2 = k 2 /t\. 

Adding and subtracting them leads to 



and 



2k 3 /h = k 1 +tl + k 2 /h, 
2t 2 = ki+tl-k 2 /h 

whose product can finally be written as 



(17) 
(18) 







fc 2 



which is a cubic equation in t\. Any positive root of 
t\ can be used to derive all real solutions of the initial 
quartic equation (negative roots would lead to imaginary 
values of ±t\. Either sign can be used to solve the factor- 
ized quartic equation. Changing the sign corresponds to 
swapping the coefficients between the first and the sec- 
ond quadratic polynomial in equation (|16[) . Descartes' 
Sign Rule [l3| can be exploited to ensure that there is 
always at least one positive root. According to the rule 
the number of sign changes of the consecutive polynomial 
coefficients is the maximal number of positive roots. Now 
one can substitute t\ by — 1\ to determine the maximal 
number of negative roots. Since k 2 is real the constant 
coefficient — k\ is negative. The leading monomial has 
also a negative coefficient. Thus there can be two or zero 



2s? - 9sis 2 + 27s 3 
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are determined. If r 2 < q 3 the cubic equation has three 
real roots which can be found by computing 



9 = arccosr/yg 3 . 
The three roots are then given by 

„ ^ ,9 s\ 
z 1 = -2Vgcos(-) - y, 

9 + 2n si 
z 1 = -2V?cos( — - — 



Z\ = —2^/qcos(- 



— 2n si 



Their first appearance goes back to Francois Vieta who 
published them in 1615. In the case of r 2 > q 3 there is 
only one real solution and defining the auxiliary variables 



(-r + Vr 5 - <7 3 ) 1/3 



and 



V = (-r- yjr 2 - q 3 ) 1 ' 3 



allows to express the real solution simply in terms of u 
and v as 

Z\ = U + V . 
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